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We consider synchronization of coupled dynamical systems when different types of interactions 
are simultaneously present. We assume that a set of dynamical systems are coupled through the 
connections of two or more distinct networks (each of which corresponds to a distinct type of 
interaction), and we refer to such a system as a dynamical hypernetwork. Applications include 
neural networks formed of both electrical gap junctions and chemical synapses, the coordinated 
motion of shoals of fishes communicating through both vision and flow sensing, and hypernetworks 
of coupled chaotic oscillators. We flrst analyze the case of a hypernetwork formed of m = 2 networks. 
We look for necessary and sufficient conditions for synchronization. We attempt at reducing the 
linear stability problem in a master stability function form, i.e., at decoupling the effects of the 
coupling functions from the structure of the networks. Unfortunately, we are unable to obtain a 
reduction in a master stability function form for the general case. However, we show that such a 
reduction is possible in three cases of interest: (i) the Laplacian matrices associated with the two 
networks commute; (ii) one of the two networks is unweighted and fully connected; (iii) one of the 
two networks is such that the coupling strength from node i to node j is a function of j but not of 
i. Furthermore, we define a class of networks such that if either one of the two coupling networks 
belongs to this class, the reduction can be obtained independently of the other network. As an 
example of interest, we study synchronization of a neural hypernetwork for which the connections 
can be either chemical synapses or electrical gap junctions. We propose a generalization of our 
stability results to the case of hypernetworks formed of m > 2 networks. 



pacs 05.45.Xt Synchronization, coupled oscillators 
pacs 05.45.Pq Chaotic systems 
pacs 89.75.-k Complex systems 
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I. INTRODUCTION 



Synchronization of coupled dynamical systems has been the subject of a considerable amount of research (see e.g., 



6l-|llj to pinning control 
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One case 



l|-[5[) with applications ranging from adaptive synchronization strategies 
of interest is that of complete synchronization that occurs when the individual systems, if appropriately coupled, 
converge on the same time-evolution. Complete synchronization can be observed in the presence of selective coupling, 
i.e., the systems are coupled through the connections of a network. A common underlying assumption is that the 
interactions among the systems are all of the same type. For this case, it has been shown that stability of the 
synchronized state depends on the details of the underlying network topology. 

In this framework, the master stability function (MSF) approach [2] to synchronization of networks of coupled 
identical dynamical systems has been widely investigated in the literature An outstanding problem is how to 

obtain a reduction of the stability problem in a MSF form when the set of coupled dynamical systems simultaneously 
interact through different networks, with each network being associated with a distinct coupling function. 

In this paper, we will focus on complete synchronization and we will retain selective coupling but we will allow for 
different types of couplings between the systems. We assume that all the connections that correspond to the same 
type of coupling form a network and the systems are connected by more than one network. This case is relevant to 
any situation where the individual units are allowed to interact through different types of coupling. For example, 
neurons in the brain are connected through both electrical gap junctions and chemical synapses, see e.g., 2^ 21 1. 
The coordinated motion of shoals of fishes depends on the sensory capabilities of each individual fish. Fishes typically 
use vision but also chemical/fiow sensing in order to localize their mates and coordinate their individual motion with 



respect to the shoal 
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23j (as in other animal species, the number of neighbors that can be simultaneously sensed 
by each fish is typically bounded and depends on the specific kind of interaction i24j). Another example is that 
of interdependent networks, such as e.g., the coupled infrastructure of power stations and internet communication 



servers 



25| . In recent years, the possibility of cascades of faults through coupled interdependent networks has been 



pointed out as a crucial aspect with respect to the assessment and design of critical infrastructures 
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In this paper, we consider that a set of identical dynamical systems Xi = F{xi{t)), i = 1^2, N , are coupled 



through the connections of m different networks, and we refer to such a system as a hypernetwork, see e.g. 



5J| . We first consider the case of to = 2 networks (a generalization to the case of to > 2 networks will be presented 



in Sec. IV). The systems are then coupled as fohows, 

JV N 

Mt) = F{x,{t)) +a^Yl ^^Ad^jit - rg)) - Gix.{t - Tg))] + ^ B,,[Hix,{t - r^)) - H{x,{t - r,.))], (1) 

i = 1,2,...,A^, where Xi{t) = (t), is the 71-dimensional state of node i, F : R" — > i?" represents 

the dynamics of each individual unit, G : R" i?" and H : i?" — ^ R" are different coupling functions, Tg and r^, are 
(possibly) different interaction delays, and are two scalar coefficients. As can be seen from ([Ij, the interactions 
between the individual units are those of two distinct networks, which are represented by the two distinct adjacency 
matrices A = {Aij} and B = {-By }. Thus Eqs. ([1]) describe a hypernetwork of coupled dynamical systems. 
An equivalent way of writing Eqs. ([1]) is the following, 

AT N 

x,{t) = F{x.{t)) + LtjG{x,{t - r,)) + fi^ J] Lf^H{x,{t ~ m)), (2) 

i ~ 1,2, N, where — Aij — Sij J2j ^ij ^^'^ ^fj = ^ij ^ ^ij X]j ^^'^ Laplacian matrices. Say and 

{A.f } the set of eigenvalues associated respectively with the two matrices and . By construction, both matrices 
and have one eigenvalue, A;^ = and = 0, with associated eigenvector [1, 1, 1]. The nN dimensional 
state space of the system in Eqs. ([2|) contains an n-dimensional synchronization manifold T, 

Xi{t) ^ X2{t) ^ ■■■ ^ XN[t). (3) 

Note that if a solution belongs to T over a time interval [t^, to + Tmax], where Tmax ~ nia-Xr^.Tft, then the solution will 
belong to I, for any time t > to + Tmax- In this case, the synchronized solutions xi{t) — X2{t) = ... = XN{t) = Xs(t) 
is characterized by the same dynamics as that of an uncoupled system, 

Xs{t) ^ Fix^it)). (4) 

The main goal of this paper is to study linear stability of the synchronous solution (|3l4p for the set of equations 
The same problem for the case that the systems are coupled through the connections of only one network, i.e., 
= in Eq. has been intensively studied in the literature, see e.g., [2, 



30l-l36l| . For this case it can be shown 
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that linear stability of the synchronous solution can be analyzed in terms of the following low-dimensional equation 

SS:{t) = DF{x,{t))Sx{t) + a^X^DHix^t - Th))Sx{t - th), (5) 

where DF (DH) represents the Jacobian matrix of the function F (H). In particular, the condition for stability is 
that the maximum Lyapunov exponents [55] associated with Eq. ([5]) are negative for k = 1, (N — 1). Eq. (5) for 
k = N yields 

Sx{t) = DF{xs{t))Sx{t), (6) 

which corresponds to the linearized equation for the evolution in the synchronization manifold ([3]). ([5]) is a system 
of n scalar differential equations as opposed to the linearized system which is described by nN scalar differential 
equations. Hence, system ([5]) is termed low-dimensional. The nice thing about this approach is that it provides 
necessary and sufficient conditions for synchronization. Similar conditions have been obtained for networks of groups 



17| , for adap tive synchronization of complex networks 



network 
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38[ , for the pinning control problem applied to a complex 
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421. In this 



40|, and for the case that slight deviations from nominal conditions are present 
paper, we attempt at obtaining a condition in terms of a low-dimensional equation for the more complex case that 
the systems are coupled through the connections of two different networks (Eq. However, as we will see, our 

proposed problem is not easy to solve in general. 

In what follows, we first consider the case that the two matrices A and i? in ([1]) are arbitrary and we show that the 
stability problem does not admit a solution in a low-dimensional form. Then we focus on three examples of interest 
for which we show that such a reduction is possible: 

• The two Laplacian matrices and commute. 

• One of the networks (either A or B) is unweighted and fully connected. 

• One of the two networks (say e.g.. A) is such that Aij = aj, i, j ~ 1, N. 

The rest of the paper is organized as follows. In Sec. II we attempt at obtaining necessary and sufficient conditions 
for stability of the synchronous solution for a hypernetwork Yet, we show that unfortunately it is not always 
possible to reduce the problem in a low-dimensional form. However, we analyze three cases of interest for which such 
a reduction is possible. Furthermore, we define a class of networks such that if one of the two coupling networks 



belongs to this class, the reduction can be obtained independently of the other network. Numerical simulations are 
shown in Sec. III. In Sec. IV we generalize our results to the case of hypernetworks formed of m > 2 networks. A 
more general class of hypernetworks that are not described by the set of equations ([IJ is discussed in Sec. V, where 
the example of a network of neurons connected by both electrical gap-junctions and chemical synapses is presented. 
Finally, the conclusions are given in Sec. VI. 

II. STABILITY ANALYSIS 

We consider stability of Eqs. ^ about the synchronous solution Linearization of Eqs. Q about ^ yields, 

N 

Sx,{t) = DF{xs{t))5x,{t) +<J^Y. LfjDG{xs{t ~ Tg))Sxj{t - Tg) 

N 

+ <y^Yl Lf,DH{xS ~ Th))5xj{t ~Th), 



(7) 



1, 2, N . The set of equations ([7]) can be rewritten in vectorial form as follows, 

5x{t) = In® DF{xs{t))5x{t) + cr^L-^ ® DG{x,{t - Tg))5x{t - Tg) 

+ a^L^ (g) DH{xs{t - Th))5x{t - th), 



(8) 



where 5x{t) — [5xi{t)'^ ,5x2{t)'^ , (5a;Ar(i)-^]-^ and the symbol ® indicates the direct product or Kronecker product. 
Now we proceed under the assumption that at least one of the two Laplacian matrices, say L^, is diagonalizable, 
i.e., = yA"*y~^, where A"^ is a diagonal matrix with the elements on the main diagonal being the eigenvalues 
A;^^, A^, A;^ and y is a matrix whose columns are the associated eigenvectors, vi,V2, vn. Then, by introducing 
the change of variable, ri{t) = V^^ (8) InSx{t), Eq. (|8]) becomes, 

r]{t) = /at ® DF{xs{t))ri{t) + cr-^A^ ® DG{xs{t - Tg))Tj{t - Tg) 

(9) 

+ cr^S (g) DH{xs{t - Th))r]{t - Th), 

where the matrix S = V^^L^V. It would be nice if the matrix S were diagonal but unfortunately there is no 
guarantee that this will be the case in general. Then we see from Eq. ^ that, different from the classical master 
stability function derivation [2], it is not possible to decouple Eq. ([9]) in blocks, each one independent of the others. 



FIG. 1: An example of two graphs with associated commuting Laplacian matrices, (a) All the links have associated weight 
equal one. (b) All the links have associated weight equal one except for the link in black having associated weight 2 and the 
links represented as dashed arrows having associated weight -1. 

A. The case that the two matrices and commute 

A special case is when the two matrices and commute. Two matrices that commute have the property of 
sharing the same set of eigenvectors, i.e., assuming that they are both independently diagonalizable, it is possible to 
write = Vk^V-^ and ^ Vk^y-^, where A^ is a diagonal matrix with the elements on the main diagonal 
being the eigenvalues of the matrix . Thus for this case the matrix S coincides with the diagonal matrix A^ as 
S = V~^VA^V^^V = A^. It follows that equation ([9]) can be decomposed in N blocks independent of each other, 

f,k{t) = DF{xs{t))m{t) + a^X^DG{xs{t - Tg))r]k{t - Tg) + a^\^DH{x,{t - Th))vk{t - th), (10) 

k = l,...,iV, where and are respectively the (complex) eigenvalues of the matrices and i^, which are 
associated with the same eigenvectors, i.e., such that L^Vk = A^Ufc and L^v^ = X^^k- Recall that the eigenvalues 
— — ^'^'^ the corresponding eigenvector is [1, 1...1]. Then for k — N, Eq. ([TO)) yields, 

f,N{t) = DF{xs{t))m{t), (11) 

which corresponds to perturbations in the direction tangent to the synchronization manifold Q and as such are 
not relevant in determining stability of the synchronous solution. Thus a necessary and sufficient condition for 
synchronization is that the Lyapunov exponents associated with Eq. IjlOp are negative for k = 1,2,..., [N — 1). 
We now introduce a parametric equation 



f,{t) = DF{xs{t))7]{t) + vDG{xs{t - Tg))r]{t - Tg) + zDH{x,{t ~ Th))v{t ~ r^), 



(12) 



FIG. 2: A hypernetwork formed of a fully connected graph (thin black arrows) and a superimposed network of 9 directed links 
(thick gray arrows). All the links (those associated with either one of the networks) have associated weight equal one. 



where y and z are two complex parameters. We associate a master stability function with Eq. (1121) . 

M{y,z), (13) 

which returns the maximum Lyapunov exponent of Eq. (jl2p as a function of the pair of complex arguments {y, z). 
Then given any hypernetwork ([2]), stability of the synchronous solution can be evaluated by checking that AA{y, z) < 0, 
for {y, z) = (ct^A^, CT^Af ), fc = 1, 2, {N - 1). Alternatively, a necessary and sufficient condition for stability of the 
synchronized evolution is that the pairs {a^X^ X^), k — 1,2, {N — 1) fall in the region of the domain of the 
master stability function Ai{y, z) for which < 0. A similar result for the case of a single network whose topology 
is allowed to evolve in time has been previously obtained in 43|. 



However, we note that the case that the two matrices and commute is quite specific and not very likely 
to occur in practical situations. An example of two graphs with associated commuting Laplacian matrices is shown 
in Fig. 1. In Sections IIB and IIC, we present two examples for which a reduction of the stability problem ([T]) in a 
low-dimensional form is possible, even if the two matrices and do not commute. 
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B. The case that one of the two networks is unweighted and fully connected 

We consider the case that one of the two networks is unweighted and fully connected. Without loss of generality 
we take this matrix to be A, 



1, for i,j = l,...,iV, j ^ i. 
0, for i = j. 



(14) 



Then = (1 — SijN), where Sij is the Kronecker delta. An example of such a hypernetwork is shown in Fig. 2. 
We consider again stability of the synchronous solution ([3|). In what follows, we obtain a master stability function 
by only diagonalizing the (N — 1) dimensional subspace of transverse perturbations without worrying about the fact 
that these may couple into the remaining direction (which is tangent to the synchronization manifold) . 
The matrix can be diagonalized as ~ VA-^V^^, where A'^ is the following diagonal matrix, 



A^ = {A^.} 



-NO •■• 
-iV ••• 



■■■-NO 

oy 



We now look at Eq. It can be shown that the matrix ^ = F ^L^V, S = {^ij}, has the form 



— (JV-l)l — (JV-1)2 

\^ Hjvi Hjv2 



2i(jv-i) 

S2(JV-1) 

-(jv-i){jv-i) 

2jv(jv-i) oy 



(15) 



In fact, the matrix L^V has a column whose elements are all zero. This is due to the properties (i) that the sum of 
the elements in each row of the matrix equals zero, and (ii) that the matrix V has a column (the same column 
where the eigenvalue of A"^ is) whose elements are all the same. It immediately follows that V^^L^V has a column 
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whose elements are all zero. Therefore, Eq. ^ can be re-expressed as, 

i^'{t) =lN-i®DF{xs{t))7j'{t) 

- a^NlN-i ® DG{x,{t - Tg))7j'{t - Tg) (16) 



i)N{t) - DF{xs{t))m{t) - <J^DH{x,{t)) J2 ^NjViit - ^h), 

3 = 1 



(17) 



where the vector rj' = [rji, TyJ, rijf_i] , and ^' is the {N — 1) dimensional square matrix. 



-11 —12 



—21 —22 



-l(JV-l) 
— 2(JV-1) 



y— (jv-i)i — {]v-i)2 • • • — (]v-i)(jv-i) y 

We note that Eq. is independent from Eq. p7)) . Hence, we term the first as the drive system and the second as 
the response system. Note that rj' corresponds to perturbations transverse to the synchronization manifold, while r]N 
corresponds to perturbations within the synchronization manifold. Thus synchronization stability is governed by Eq. 
P^ . which does not involve tjn. We diagonalize the matrix S', obtaining (TV — 1) blocks of the form. 



Cfe(t) = DF{x,it))Ckit) - <J^NDG{Xs{t - Tg))Ck{t - Tg) 

+ a^VkDH{xs{t - Th))Ck{t - Th), 



(18) 



k = 1, {N — 1), where (I'l, V2, vn-i) are the eigenvalues of the matrix S'. Note that the eigenvalues of the matrix 
S' are the same as those of the matrix L^, except for the one eigenvalue that is equal to 0. 

If the {N — 1) maximum Lyapunov exponents associated with the drive system (|18p are all negative, then for large 
enough t, Cfe(i) — > 0, fc = 1, [N — 1). If this happens, then Eq. (flT)) yields for large enough t, 



(19) 
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which corresponds to the hnearized equation in the direction tangent to the synchronization manifold. 

Thus we can introduce the parametric equation ()12|) in the pair {y, z), with the parameter z being possibly complex 
and a master stability function (jl3p which returns the maximum Lyapunov exponent of Eq. (|12p as a function of 
the parameters y and z. For a given hypernetwork (|2ll4p . a necessary and sufficient condition for stability of the 
synchronous solution (jS]), is that y = —a^N and z — u^Vk, k — 1,2,..., {N — 1) belong to the region of the domain 
of the master stability function p^ . for which M{y, z) < 0. 

This formulation allows to decouple the effects of the dynamical function F and the coupling functions G and 
H, from those of the network matrices and . In particular, for any given triplet of functions F,G, and H, 
the matrix B determines the parameters 1^1,1^2, i^N-i, and if the master stability function A4{y, z) is negative for 
y ~ ~a^N and z = cr^t^i, <j^U2, cr^j^(jv-i)i then the synchronization manifold is stable. An interesting thing about 
our derivation (|18p is that we have been able to obtain a reduction of the stability problem ([T]) in a low-dimensional 
form though the two matrices and do not necessarily commute. 



Here we consider the case that the coupling strength from node j to node i is only a function of the source node j 
and not of the destination node i, that is 



An example of such a network is shown on the left-hand-side of Fig. [3l where the width of each link j — > i represents 
the strength of the associated coupling Aij. The network on the right-hand-side of Fig. [3] is an outward star graphs 
corresponding to setting Oj — for j = 1, {N — 1) and a^r 7^ in Eq. ((20|) . Under assumption (|20|) . Eqs. ^ 
become, 



C. 



The case that Aij — aj 




i,j = l,...,N. 



(20) 



N 




(21) 



N 



+<j^ [H{x,{t - T,,)) - H{x,{t - Th)) 
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1, 2, .., iV, which can be recast in the form of Eq. ([2]), with the matrix = having the form. 



^ a\ — a a2 



a\ a2 — a 



y 0.1 a2 



a(jv-i) cln 
O(N-i) ~ <i On 



(22) 



where a = Y^^=i^j- The matrix in ([22|) has the property that it has one eigenvalue = with associated 



eigenvector [1, 1, 1] and the remaining [N — 1) eigenvalues are = A 
diagonahzed as = VK'^V^^ with equal, 



-a. Moreover can be 



... n\ 



-a ... 

-a ... 

... -a 

oy 



It is easy to see that the matrix S is in the form (fT5|) . with the entries in the A^-colunin being all equal to zero. This 
allows to decouple the set of linearized equations into a drive subsystem and a response subsystem, with the response 
subsystem corresponding to perturbations tangent to the synchronization manifold ^ and the drive subsystem 
corresponding to perturbations transverse to the synchronization manifold. 

Then following the analysis in Sec. IIB, it can be shown that a necessary and sufficient condition for stability of 
the synchronous solution for the hypernetwork (|2ip is that the maximum Lyapunov exponent of the low-dimensional 
equation, 



(23) 



is negative for k ~ 1, {N — 1), where Af , A^, A^_-^-j are the eigenvalues of the matrix , excluding the one 
eigenvalue A^ = 0. It is then possible to associate Eq. (1^^ with the parametric equation (fT^ and the master stability 
function (jl3p , which returns the maximum Lyapunov exponent of Eq. (jl2p as a function of the pair of parameters 



{y,z), with the parameter y — —a^a and the (possibly complex) parameter z = cr^Af , fi^Af , i7^A^_-^^. Again 
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FIG. 3: On the left-hand-side, a. N = 5-node network belonging to class C, i.e., such that the entries of the associated adjacency 
matrix A = {^ij} satisfy Aij = aj (the condition discussed in Sec. IIC). The width of each link j — ^ i represents the strength 
of the associated coupling Aij. The network on the right-hand-side is an outward star graph, corresponding to satisfying (|20p 
with aj = 0, j = 1, {N - 1). 

we note that we have been able to obtain a reduction of the stabihty problem ([7]) in a master stability function form 
though the two matrices and do not necessarily commute. 

We wish to emphasize that the case in Sec. IIB (fully connected network) can be seen as a subcase of that in Sec. 
IIC {Aij = aj). In fact, if we assume aj = a, j = 1, ...,N in (I20|) . then the Laplacian matrix L"^ = {Lfj} in ([22l) is 
such that Lfj — a(l — SijN), i.e., it coincides with the matrix L'^ considered in Sec. IIB up to a multiplicative factor 
a. 

D. Necessary conditions on the matrix A 

We observe here that there is a substantial difference between the conditions on the adjacency matrices A and B 
(the Laplacian matrices and L^) discussed in Sec. IIA and those discussed in Sees. IIB and IIC. First consider 
the case presented in Sec. IIA, that the two Laplacian matrices and commute; then, if one of the two matrices 
changes, there is no guarantee that the condition would still hold. On the other hand, the conditions discussed in 
Sec. IIB and IIC refer essentially to one of the two matrices, allowing the other one to be freely chosen. 

In sections IIB and IIC, we have found sufficient conditions on one of the two adjacency matrices, say A, that if 
satisfied, allow a reduction of the stability problem in a low dimensional form, irrespective of the other adjacency 
matrix, say B. In this section, we are interested in finding necessary condition for this to happen. We consider the 
set of Eqs. ((T|) and we define the class C of all the networks A that satisfy the property of allowing a reduction of 
the stability problem in a low dimensional form, irrespective of the other network B. In what follows, we show that 
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a network in C is such that the entries of the associated adjacency matrix A = {Aij} satisfy Aij = aj, i.e., the same 
condition discussed in Sec. IIC. 

Hereafter, we seek to find the conditions for an adjacency matrix A (a Laplacian matrix L^) to be in C. Based on 
our previous discussion in Sees. IIB and IIC, we see that the properties that the matrix has to satisfy are the 
following: 

(A) is diagonalizable. 

(B) The sums of the elements in the rows of the matrix are equal zero. This also implies that the matrix 
has one eigenvalue equal zero, with associated eigenvector [1,1,...,!]. 

(C) The remaining {N — 1) eigenvalues are all the same. 

If the three properties above are satisfied, the matrix can always be written as follows, 

= WPW-\ (24) 

where the matrix P is a diagonal matrix with all the entries on the main diagonal being equal to the same value, say 
p, except one entry (which, without loss of generality, we assume to be the one in the rightmost column) that is equal 
to zero. The matrix W is any invertible matrix with the rightmost column being equal to the vector [1, 1, 1]. We 
note that the matrix P can be rewritten as P = p{In — 1^)^ where In is the identity matrix and is a diagonal 
matrix with all the entries on the main diagonal being equal to zero except the one in the rightmost column being 
equal to one. It follows that 

L^=p{I-WI*nW-^). (25) 

It is easy to see that the matrix WI'^W~^ is by construction such that the entries in each one of its columns are the 
same. Hence, the corresponding adjacency matrices A have to be in the form Aij — aj, discussed in Sec. IIC. 

We conclude that if we are given a specific adjacency matrix B (a specific Laplacian matrix L^), there are two 
possible choices of the adjacency matrix A (the Laplacian matrix L^) for which the stability problem can be reduced 
in a low-dimensional form: (i) commutes with , and (ii) A belongs to C, i.e., its entries are such that Aij — aj. 
Note that condition (ii) is independent of the choice of the matrix . 
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III. EXAMPLES 



A. Example 1: Coordinated motion of swarms of particles. 



Swarms of birds, hordes of insects, shoals of fishes, and colonies of ants have been modeled as systems of interacting 
self-propelled particles 



23 



44 



45|. Here we consider a simple model of N particles moving along a fixed direction, 
say y, through a resistent fluid. The position (velocity) of particle i along the y direction is labeled as yi{t) {vi{t)), 
i — 1, ...,N. We consider the following equations of motion. 



m (26a) 
v^it) =(a - /3<(t)2)<W +J2m,{y,{t) - y,{t)) + ^ m,c,,{t){v,{t) - v,{t)), (26b) 

3 3 

1 = 1, ...,N. The first term on the right hand-side of Eq. (|26b ') represents propulsion/friction of particle i, w[(i) is 
the relative velocity along y with respect to the resistent fluid of particle i, vl{t) — {vi{t) — Vf{t)) and Vf{t) is the 
velocity of the resistent fluid, which we model as an external input and we assume to be uniform in space. The second 
term on the right hand-side of Eq. (|26b ) represents attraction from particle j on particle i. The third term on the 
right hand-side of Eq. (|26b ) models a relative velocity adjustment between particles, rrij > is the mass of particle 
j — 1, N, a, /3 > 0, Cij(t) measures the strength of the interaction from particle j on particle i, which we set to be 
a function of the physical distance between particles i and j. 



c.,(^) = [4 + (%W-y.Wf]^ (27) 

where dij is the distance between particles i and j in the plane orthogonal to the y direction and the exponent e 
determines the strength of the interaction as a function of the distance. An analogous model for particles that are 
allowed to move in the three-dimensional space has been considered in |46| . 

We note that the system of equations (1^5)) admits a synchronous solution yi{t) — y2{t) = ... = yN{t) = ys{t), 
vi{t) = W2(t) = ... = Ujv(i) = Vs{t), obeying 



vsit) =[a - l3Mt) - vf{m{vs{t) - vf{t)), 



(28a) 
(28b) 
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where again is an external input. The synchronous solution corresponds to a configuration in which all the 

positions and the velocities of the particles along the y-direction are the same. We are interested in studying stability 
of this solution. In order to do that, we linearize Eq. (l26l) about (pS]) . 



5yS) =Sv,{t), (29a) 
6v,it) =[a - 3/?K(i) - it))^]6v,{t) + ^ (.Jy, (t) - 6y,it)) + ^m,K,)''('5i^j (0 - Sv^it)), (29b) 



i — 1,...,N. Equations (j29p can be rewritten in matrix form, 



1 

[a ~ 3(3(vAt) - vf {t))^] 



where 



^ ^A,,[fe,(t)-<5x.(i)]+ ° I J2B,,[5x,it) -Sx.it)], 
i \ 1 / J 



(30) 



Sxi{t) = 




(31) 



and Aij — nij, Bij — mj{dij)'^'^, i,j — 1, N. It is easy to see that the matrix A — {Aij} belongs to class C. Hence, 
following Sec. IIC, the stability problem can be reduced in a low-dimensional form analogous to Eq. 



ykit) = I 1 9k{t), 

[a-3f3{vs{t)-vf{t)f] + X^ 



(32) 



fc = 1, (N— 1). Note that a — > 0- Thus a necessary and sufficient condition for the synchronous solution 

to be stable is that < (ws(t) — {t))"^ >t> {a + A^)/(3/3), fc = 1, {N — 1), where again we have used the symbol 
< ... >i to indicate the time- average. 
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FIG. 4: Sign of the master stability function M{y,z) for a network of Rossler systems ((33}, G{x(t)) — [xi{t),0,0]'^ and 
H(x{t)) = [Q,X2{t),Qf. The gray (white) area indicates a negative (positive) maximum Lyapunov exponent. 



B. Example 2: Synchronized Chaotic Motion. 



In what follows, we consider a hypernetwork that allows a chaotic synchronous evolution We choose rt = 3, 



F{x) 



-~X2 - ^3 

xi + 0.2x2 
0.2 + (xi - 7)x3 



(33) 



is the equation of the chaotic Rossler oscillator, G{x{t)) — [a;i(i), 0, O]"'", and H{x{t)) — [0, a;2(i), 0]"^, Tg = Th = 0. 
Stability of the synchronous solution for networks of coupled Rossler oscillators coupled via different coupling functions 

n 

has been widely investigated in the literature, see e.g., (IS^. While it is known that this problem allows a low- 
dimensional reduction, the case of dynamical hypernetworks has not been considered. In what follows, we show that 
under specific conditions, a low-dimensional analysis can still be applied and based on this approach, we derive new 
conditions for the stability of the synchronous solution. We numerically compute the master stability function M{y, z) 
associated with Eq. (fT2)) for the case that y and z are real numbers. 

Fig. 4 shows the results of our computations with the gray (white) area indicating a negative (positive) master 
stability function. We wish to emphasize that once the master stability function has been computed for a given triplet 
F, G, and H (as shown in Fig. 4), we are able to predict stability of the synchronous solution for any hypernetwork 
([2]), corresponding to either one of the three cases presented in Sections IIA, IIB, and IIC. 
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LU 5 



(a) 

1 




1 
1 

1 






FIG. 5: (b) shows the intersections between the right profile of the stability area of Fig. 3 and the lines y — —4, y = —4.5, 
and y = —5.5. (a) and (c) show the results of numerical simulations for which we have have integrated Eqs. and (|33p with 
G{x{t)) — [a;i (t), 0, 0]"^ and H{x{t)) — [0, a::2(t), 0]"^ for a long time and recorded the average synchronization error E. 



We define the average synchronization error E, 



N 



(34) 



1=1 £=1 



where Xie{t) — N^^J2^=i^ie{'t): Pi = < {xsi— < Xsi >)^ >^/^, < ... >t indicates a time average and 
Xs = {xis,X2s,X3s)^ denotes the dynamics of an uncoupled system (i.e., using dynamics from Eq. (HJ). 

We consider the hypernetwork shown in Fig. 2. We assume that the matrix A is associated with the fully connected 
network (whose connections are represented as thin black arrows in the figure) and that the matrix B is associated 
with the superimposed graph (in gray in the figure), i.e., the entries of the matrix B are Bij = 1 if there is a direct 
arrow from node j to node i in the figure, Bij = otherwise. Then we have that the eigenvalues of the matrix 
are {—3, —2.618, —2, —1, —0.382, 0}; note that they are all real and less or equal zero. In general, in order to verify 
stability, it is necessary to check that all the pairs (y, z) = (o"^A^, cr^ X^), k = 1,2, {N — 1) follow into the domain 
of the master stability function for which A4{y, z) < 0. This can be done for example by superimposing the {N — 1) 
points corresponding to all the pairs {a^X^ X^) to Fig. 4; if all the points fall into the grey area, this ensures 



stability (sufficient condition for synchronization) and if only one of the points fall into the white area, this corresponds 
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to instability (necessary condition for synchronization). However, for the network of Fig. 2 and the master stabihty 
function of Fig. 4, we note that for any fixed value oi y = —a^N, the condition for stability is that 

a^Af</c, i = l,...,(iV-l), (35) 

where the parameter k is the abscissa of the intersection of the y = —a^N line with the right profile of the gray 
area shown in Fig. 4. Note that n can be either positive or negative. We define X^^ix — "^sx (Af , Af , A^ -^-j) and 
A^j„ — min (Af , Af , A^ -^-j). Then for this case, stability of the synchronous solution can be assessed by testing 
the following simple condition, 

<^^^?nax <l^, if K < 0, 
<^^^Ln < if K > 0. 

In Fig. 5 we consider the following three cases: cr^ = 4.5/6, = 5.5/6, and <t^ ~ 2/3 (corresponding respectively 
to, y = —4.5, y — —5.5, and y — —4). As can be seen from Fig. 5(b), for the first two cases k < 0, while for the latter 
K > 0. We integrate Eqs. @ and §^ with G{x{t)) = [xi(t), 0, 0]"^ and H{x{t)) = [0,X2{t),0f for a long time and 
record the average synchronization error E. As can be seen from Fig. 5(a) (5(c)), E approaches zero iff <J^\^ax 
when CT^ — 4.5/6 and — 5.5/6 {<J^ ^min < when ct^ — 2/3), thus confirming the master stability functions 
predictions. 

IV. GENERALIZATION TO m NETWORKS 

In this section we consider synchronization of a hypernetwork formed of m > 2 distinct networks. For this case, we 
rewrite Eq. ^ as follows, 

m N 

xdt) = F{x,{t)) +Y.<r^Y. L%GHx,{t - r'^)), (37) 

k=l 3 = 1 

« = 1, 2, N , where : i?" — ?> i?" is the coupling function associated with the connections of network fc, — {ijj} 
is the Laplacian matrix associated with network fc, is a scalar measuring the overall coupling strength for network 
k, k = \, ...,m. In what follows, we will generalize the main results of Sec. II to this more general case (Eq. ([57)) ). 
The delays r*^ may be possibly different, i.e., ^ Tj, i,j — 1, ...,m, i ^ j. The nN dimensional state space of the 
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system described by Eqs. (p7)) contains the n-dimensional synchronization manifold T, defined by Eq. ([3]). Note that 
if a solution belongs to X over a time interval [io, + Tmax], where Tmax = max^ r*, then the solution will belong to I, 
for any time t > to +Tmax- In this case, the synchronized solutions xi{t) = X2{t) = ... = x^it) — Xs{t) is characterized 
by the same dynamics as that of an uncoupled system Q. In what follows, we are interested in evaluating stability 
of the synchronization manifold I. 

As a first case, we consider that the matrices {L'^}, k — 1, ...,m all commute with each other, i.e., they all share 
the same set of linearly independent eigenvectors. Then, similar to Sec. IIA, it can be shown that stability of the 
synchronous solution can be reduced in the following low-dimensional form, 



I = 1,...,A^, where {Af} is the set of (complex) eigenvalues of the matrices {L''}, which are associated with the 
same eigenvectors, i.e., such that L'^vi — Xfvi, k = 1,...,to and I = 1,...,N. Recall that for any k — l,...,m, the 
eigenvalue = 0, and the corresponding eigenvector is [1, 1...1]. Hence, for k — N, Eq. (1551) yields Eq. dTTI) which 
corresponds to perturbations in the direction tangent to the synchronization manifold ([3]) and as such are not relevant 
in determining stability of the synchronous solution. Thus a necessary and sufficient condition for synchronization is 
that the Lyapunov exponents associated with Eq. ([55)) are negative for k = 1,2, (N — 1). It is then possible to 
associate the following master stability function with Eq. (|38p 



which returns the maximum Lyapunov exponent of the system psp for — a^X^ . A necessary and sufficient condition 
for stability is that y™) < for / 1,...,(7V- 1). 

We now attempt to generalize the result of Sec. IIC to a hypernetwork formed of m networks. We assume that the 
first (to — 1) networks, fc = 1, (to — 1), belong to C, while the remaining network, fc = m, is arbitrary. Under these 



m 




(38) 



M{y\y^....yn. 



(39) 
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assumptions the first {m — 1) Laplacian networks are in the following form: 



ai — a' a2 



a-i an — a 



k 
0-1 



k 
0,2 



k 

O(N-l) 



'■(N-1) 



k 



*(iV-l) 



(40) 



where a 



1). Note that two matrices in C, i.e., having the form (PO]) . do not necessarily 



commute. Each matrix L*' in (|40| has the property that it has one eigenvalue = with associated eigenvector 



[1, 1, 1] and the remaining [N — 1) eigenvalues are \\ = Af = ... = A^^^ ^^j ~ —a'^, k ~ 1, {m — 1). 

The eigenvectors of any of these matrices can be used as a new basis, say we choose k — 1, — VA^V^^. 
Then it is easy to see that all the matrices V^^A'^V, for k = 2,...,m, are in the form ((T5|) . It follows (similarly to 
Sec. lie) that we can decouple the set of linearized equations in a drive subsystem and a response subsystem, with 
the response subsystem corresponding to perturbations tangent to the synchronization manifold ([3]) and the drive 
subsystem corresponding to perturbations transverse to the synchronization manifold. Moreover, it can be shown 
that a necessary and sufficient condition for stability of the synchronous solution for the hypernetwork (1211) is that 
the maximum Lyapunov exponent of the low-dimensional equation. 



k=l 



is negative for I = 1, (A'^ — 1), where A™, A™, are the eigenvalues of the matrix L"*, excluding the one 
eigenvalue A^ = 0. A necessary and sufficient condition for stability is that the master stability function ([M)) is 
negative for I = 1, (TV - 1), where y'' = -a'^d'', fc = 1, (to - 1) and = cr"A[", I = 1, {N - 1). 

Finally, we consider the more general case that m' < m networks of the hypernetwork (j37p belong to C and the 
remaining (m — m') Laplacian matrices commute with each other. Without loss of generality, we assume that the 
first to' networks in p7p are in C, A: = 1, ...,m', and that the remaining (to — to') Laplacian matrices L'' commute 
with each other, k = (m' + 1), m. We observe that a reduction of the synchronization stability problem in a low 
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dimensional form is possible, 

m' m 

f,i{t)^DF{x.it))m{t)-Y,^''^''DG\xs{t-T'')Mt~T'')+ (^''>HDG\xS-T''))r]i{t-T''), (42) 

k=l fe=(m'+l) 



I = 1,...,N, where Af , A2 , A^^-^^ are the eigenvalues of the matrix L , k = (to' + 1),...,to, which are associated 
with the same eigenvectors, i.e., such that L'^vi = Xfvi, k = [m' + l),...,m and stability of the low-dimensional 
equation can be associated with the master stability function ((39)) . where 

. —cF^a^, k — 1, m', 

y' = { (43) 

cr'^Af , k ^ (to' + 1), ...,m, 

I = 1, (iV— 1). The eigenvalue A™ = ... = A™ = 0, with associated eigenvector [1, 1]-^, represents perturbations 
tangent to the synchronization manifold and as such is not relevant in determining stability of the synchronous solution. 

V. STABILITY ANALYSIS FOR A MORE GENERAL CLASS OF HYPERNETWORKS 

In this section, we consider hypernetworks of coupled systems, which cannot be cast into the specific form of Eqs. 
P]). We will show that under appropriate conditions, the master stability reduction studied in Sec. II can be extended 
to study synchronization for this more general class of hypernetworks. In particular, we focus on synchronization of 
neuronal networks. Global synchronization of large areas of the brain is usually associated with the onset of a 
pathological condition, such as Parkinson's disease or epilepsy 47[. 



We study a hypernetwork of neurons coupled through both chemical synapses and electrical gap junctions. Such 
neuronal networks of different types connecting the same set of neurons have recently been explicitly discussed in the 
context of the C. Elegans nervous system, which has both a gap junctional network and a chemical synaptic network 



48 



49|. Following 
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50( 1 , a neuronal hypernetwork with these characteristics can be described by the following 



system of differential equations. 



~ '^^2;i(i)) ^ AijS,y(i)cr ^^^^ 
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where the n-dimensional vector Xi{t) = [xl{t),xf{t),...,x'2{t)] is the state of neuron i, with the first variable xl{t) 
representing its membrane potential, F : R" — >■ R" defines the dynamics of an uncoupled neuron, the coupling matrix 
A = {Aij} specifies the connection topology of the network of chemical synapses j ^ i, while the coupling matrix 
B = {Bij} specifies the connection topology of the network formed of electrical gap junctions j ^ i, kf = Aij, 
fcf = Bij, and are two scalar coefficients, Ej is the synaptic reverse potential of neuron j. Note that the 
matrix A (B) is assumed to be asymmetrical (symmetrical). The n-matrix 



1 ••• 
•■■ 





specifies the form of the coupling, indicating that neurons are coupled through their membrane potentials, the n- vector 
= [1,0, ...,0]"^ has a similar function, that is, selecting the first state variable x] of the state-vector x.^; F = <;?-^. 
The dynamical variables Sij (t) represent how strongly cell j is connected to cell i and obey the following differential 



equation 



,{t) = -cis,,{t) + C2(l - s,,{t))S{^'^x,{t - r)). 



(45) 



i,j = 1,...,N, where t is the interaction delay associated with synaptic coupling (due to axonal conduction and 
synaptic processes), ci, C2 > are two scalar coefficients, 5 : i? — > i? is a sigmoidal function, which we set. 



S{c;^Xj{t-T)) = 1 + tanh{{<,'^Xj{t ~ t) - vth)/vsi). 



(46) 



where represents the slope of the function S when its argument is small and Vth is the firing threshold. As can be 
seen from (j44[) . the individual neurons may simultaneously interact through two distinct networks, i.e., the network 
A formed of chemical synapses and the network B formed of electrical gap junctions. 
The condition for the set of equations ([^^ to admit a synchronous solution 



Xi{t) =X2{t) = ... = XN{t) = Xs{t), 
Sll{t) =Si2(i) = ... = SNN{t) = Ss{t), 



(47a) 
(47b) 
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is that _Ei = i?2 = = En = Eg. If this condition is satisfied, the synchronous solution Xs{t) obeys, 

Xs{t) ^F{xs{t)) + a^{Es - c^'^xs{t))ss{t)c;, (48a) 
s,{t) = - cis,{t) + C2(I - Ss{t))[l + tanh((?^x,(t - r) - Vth)/v,i)]. (48b) 

Note that differently from the case considered in Sees. I, II, and III, the synchronous solution (l47l) does not obey the 
same equation as that of an isolated system. Our goal in this section is to study stability of the synchronous solution 
(l47l) for the hypernetwork (|44)) . In order to do that, we linearize the set of equations (l44l) about (|47|) . obtaining 

N N 

Sx^it) ^DF{xs{t)) - <j^Tss{t)]Sx^{t) + o^^{E, - ^'^ Xs{t))Y,A[^5s,,{t) + ^ B[jV[5xj{t) - 5x,{t)], (49a) 
5s,,(t) = - ciSs,,it) - C2S{<.^x,{t - r))5s„(0 + C2(l - s,[i^)DS{<.^ x,(t - t))<{^ 6x,{t - r), (49b) 

where the matrices A! ~ {^ij} and B' = {B[^ are such that A!^ — (kf) ^ Aij and = (kf) ^ Bij and we have used 
the properties that A'^^ = 1 and B[.^ — 1. We introduce the perturbation 5ai{t) — J^j * — ^, ■■■,N. 

By multiplying (j49b ) by A^^- and summing over j, we can rewrite (j49p as 

N 

6x,{t) ^DF{x,{t)) - a^rss{t)]6x,{t) + a\{E,s - x,{t))Ss,{t) + a"" ^ Bl^r[dx,{t) - 6x^{t)], (50a) 

Shit) = - ci<5s,(t) - C2Si<;^x,{t - T))6s,{t) + C2(l - s,it))DSi^'^Xsit - r))^^ ^ 4/a;,(t - r), (50b) 

j 

We can now introduce the (n + l)-vectors 5xi{t) = [5a;i(i)^, (5s,;(t)]"^, i = 1, and the {N{n + l))-vector (5i;(t) = 
[5xi{tY' ,5x2{tY' , 5iAr(t)"^]^. Then, following Sec. II, we can rewrite the set of equations (|49l) in vectorial form as 
follows, 

5i{t) =In® [DFi{xs{t),Xs{t - r), Ssit))]Sx{t) + A' (g) DF^ixsit - r), Ss{t)))5x{t - r) + L^' ® f5x{t), (51) 
where the Laplacian matrix L^' = {i^'^^ j = {B[^ — 5ij} and the (n + l)-square matrices 

DFl{Xs{t),X,(t-T),Ss{t)) = 



DF{Xs{t)) - cr^rs,(t) +o-^?(i;. - c;'^Xs{t)) 
-ci -C2S'(<?^a:,(t-r)) 



(52) 
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DF2{Xs{t~T),Ss{t)) 





i^C2il-Ss{t))DSic'^Xs{t-T)) 



(53) 



r 





(54) 



As can be seen, the structure of the hnearized equations ([51]) is quite similar to that of Eq. ([8]) in Sec. II. The main 
difference with Eq. ([5]) is that in the case above, one of the two couphng matrices, namely A' , is not a Laplacian 
matrix, as the entries along each row of the matrix A' sum to one and not to zero. We now wonder whether the 
stability problem (j5ip can be reduced in a low-dimensional form. As for the case of Eq. ([5]), the main difficulty is 
that in general it is impossible to decouple Eq. ([5T|) in N independent blocks. One possibility, which we do not 
give further consideration in what follows, is that the two matrices A' and L^' commute. Another possibility is that 
the matrix L^' belongs to class C. If this is the case, then the matrix L^' can be diagonalized as in Eq. (|24|) . i.e., 
L^' = W{I^ — In)W~^, where the matrix W is an invertible matrix with the rightmost column being equal to the 
vector [1, 1, 1] and is a diagonal matrix with all the entries on the main diagonal being equal to zero except the 
one in the rightmost column being equal to one (see Sec. IID). Under these assumptions, the matrix S = W^^A'W 
has the form 



—11 



— (JV-l)l — (JV-1)2 



v 



Hi(jv-i) 

S2(JV-1) 

-(jv-i)(]v-i) 

— jv(jv-i) ly 



(55) 



from which we see that similarly to Sec. IIC, the linearized problem (j5ip can be decoupled into a drive subsystem and 
a response subsystem, with the response subsystem corresponding to perturbations tangent to the synchronization 
manifold (1471) and the drive subsystem corresponding to perturbations transverse to the synchronization manifold. 

5l| and in the posterior part of the putamen ^52)] , small 



It is known from the literature that in the visual cortex 
groups of neurons are likely to form dense and uniform clusters of gap-junctions. Hence, assuming that the network 
L^' is of class C can be appropriate to model such agglomerates of neurons or small subsets of them. Therefore, as an 
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example, we consider a small group of N neurons connected by a dense L^' network of gap junctions, with L^' G C. 
Under these assumptions, by diagonalizing the {N — l)-dimensional subspace of transverse perturbations (see Sec. 
II), the high-dimensional problem (j49p can be reduced into the low-dimensional form, 

m - [DF,ixsit),xAt - r), - a^f]^(i) + X^'.DF^ixsit - r), s,{t)))^{t - r), (56) 

fc = 1, 2, (N - 1), where {A-^^}, fc = 1, N, are the eigenvalues of the matrix A' . By construction, the matrix A' 
has one eigenvalue, A"^^ — 1, with associated eigenvector [1, 1, 1]. This eigenvector represents perturbations that 
are tangent to the synchronous solution, hence it is not relevant in determining stability of the synchronous solution 
(gZD. 

In the more general case in which L^' does not belong to class C and the two matrices L^' and A' do not commute, 
stability of the synchronous solution results in a much more complex problem, for which (I49p cannot be reduced in a 
low dimensional form and we expect a higher degree of complexity. The study of this case is beyond the scope of this 
paper. 

We run numerical simulations in which each individual system is described by the FitzHugh-Nagumo model, n — 2, 



F{x) = 



10[xi(a;i - 0.1)(1 - a;i) - + 0.2] 
Xi — 0.5a;2 



(57) 



and we set vth = 0.3, v^i = 10~^, ci = C2 = 10, Eg = 1, (t^ = 1, t = 1. In Fig. 6(a) we plot the time evolution 
of the synchronous evolution, obtained by integrating Eq. (|48p for this particular choice of the function F in ([57)1 
and of the parameters. We further set — 0.9 and calculate the maximum Lyapunov exponent associated with 
the low-dimensional system (|56p as a function of the parameter A'^'. This corresponds to a master stability function 
(MSF), which is plotted in Fig. 6(b) for the case that its argument is real. As can be seen from Fig. 6(b), the MSF 
curve crosses the 0-ordinate line at two distinct values of the abscissa, which we found to be approximately equal to 
—0.74 and 1.1 (in the figure, the 0-ordinate and the 1-abscissa lines are plotted as dashed lines). Thus a necessary and 
sufficient condition for stability of the synchronous solution is that —0.74 < < 1.1, i = 1, (iV— 1). If we assume 
A'lj > 0, we have by the Perron-Frobenius theorem that \X^[\ < 1, i — 1,...,N, where 1 is the Perron-Frobenius 
eigenvalue of the matrix A' , and the necessary and sufficient condition for stability reduces to — 0.74 < A"^^j„, where 

A mm = ™ini=l,...,(iV-l) i- 
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FIG. 6: (a) Time evolution of the synclironized solution for the system (|44p . obtained by numerically integrating Eqs. (|48l57p 
with Vth = 0.3, Vai = 10~^, ci = C2 = 10, Ei = E2 — ... = En = Eg = 1, — 1, t = 1. xl(t) is plotted in black and 
Ss{t) is plotted in gray, (b) Plot of the master stability function corresponding to the low-dimensional system (|56p versus the 
parameter for the case that is real. The parameters are the same as in (a), = 0.9. 



We finally run simulations of the full nonlinear hypernetwork described by Eqs. (|44l45l57p . We set the initial 
conditions for x} and x^, i = 1,...,N and for Sij, i,j ~ 1,...,N to be random numbers drawn from a uniform 
distribution in the range (0, 0.2). We consider that the network of chemical synapses is the network of = 6 nodes 
and 9 directed links represented in gray in Fig. 2, i.e., the entries of the matrix A are Aij = 1 if there is a gray direct 
arrow from node j to node i in the figure, Aij — otherwise. The spectrum of the corresponding matrix A' is real 
and X^'min = — V2/2 > —0.74. We set the network of chemical synapses to be such that Bij = bj ~ j, i,j ~ 1, 6 
(note that the particular choice of the values of bj, j — 1, ...,N affects neither the spectrum of the matrix L^' nor 
the low-dimensional equation (1561) ). We evolve the hypernetwork (I44I45I57P from t = to t — 500. We monitor 
the quantity E{t), defined in Eq. As expected, we observe that after a transient, E{t) 0. We repeat the 

same experiment for the case that Aij = 1 if |« — j| = 1 and Aij — otherwise. For this case, the spectrum of the 
corresponding matrix A' is real but X^'^in = — 1 < —0.74, thus predicting that the synchronous solution is unstable. 
This is confirmed by our numerical experiments, showing that, when the full nonlinear system (|44l45l57p is integrated 
from initial conditions that are close to the synchronous state (US]), E{t) does not converge to 0. 



VI. CONCLUSION AND DISCUSSION 



In this paper we have studied synchronization of coupled dynamical systems when different types of interactions 
are simultaneously present. Our study applies to any situation where the individual units interact through different 
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coupling mechanisms. For example, neurons in the brain are known to be connected through both electrical gap 



junctions and chemical synapses, 
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48|,|49|. Also, our study encompasses a situation where different coupling 
functions correspond to different interaction delays. 

In our formulation, a set of identical dynamical systems are coupled through the connections of two or more distinct 
networks (each of which corresponds to a distinct coupling function) and we refer to such a system as a dynamical 
hypernetwork. We first focus on the case of a hypernetwork formed of to = 2 networks and we seek to obtain necessary 
and sufficient conditions for synchronization. In Sec. II we try to reduce the stability problem in a master stability 
function form. Though a solution in this form seems to be not available in general, we show that such a reduction is 
possible in three cases of interest: (i) the Laplacian matrices associated with the two networks commute; (ii) one of 
the two networks is unweighted and fully connected; (iii) one of the two networks is such that the coupling strength 
from node j to node i is a hmction of j but not of i, with case (ii) being a subcase of (iii). We introduce a unique 
master stability function that determines stability for all three cases. Also, we define the class C of networks for which 
the reduction is always possible, independent of the structure of the other network. 

We note that in many situations, such as, e.g., in biological networks, different types of interactions are typically 
present, but the couplings may vary in time due to changing environmental conditions, making satisfaction of either 
one of conditions (i), (ii), or (iii) difficult. On the one hand, this highlights a limitation of the master stability function 

n 

approach that does not seem to be applicable to situations of arbitrary complexity (see also e.g., XB)- the other 
hand, it poses the fascinating challenge of defining alternative tools to addressing stability for the case of arbitrary 
hypernetworks. We also point out here that we cannot exclude the existence of other conditions to be satisfied 
simultaneously by both matrices A and B (e.g., for either the hypernetwork ([Ij or (HU) that allow a reduction of the 
stability problem in a low-dimensional form. 

In Sec. IV we have proposed a generalization of our stability results to hypernetworks formed of m networks. In 
Sec. V we have shown the possibility of generalizing our approach to hypernetworks of coupled systems, which cannot 
be cast into the specific form of Eqs. ([iJ. As an example of interest, we have studied synchronization of a neural 
hypernetworks for which the connections can be either chemical synapses or electrical gap junctions. The results of 
this paper could be also easily extended to study synchronization of dynamical hypernetworks of coupled discrete-time 
systems. 
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Appendix: The special case of hypernetworks of N — 2 nodes 

In this appendix we show that for hypernetworks of A*" = 2 nodes, the stabihty problem can always be reduced in a 
low-dimensional form. We start by considering that N is an arbitrary number and that the hypernetwork is formed 
of m = 2 networks (Eq. ([T|)). The generalization to the case of m > 2 networks is straightforward. 

We look at Eq. ([2]). In general, a case of interest is that one of the two Laplacian matrices, say L^, can be rewritten 

as, 

= kiL^^ + k2L^^, (58) 

where the matrix L"^^ belongs to C (i.e., it is in the form of the matrix (f22|) ) and the matrix L^'^ commute with , 
that is ^ VA'^V-^ and = VA^V-\ where A-^ and A-^ are diagonal matrices. Under the condition (|58p . Eq. 
dS]) can be rewritten as, 

Sx{t) ^/at (g) DF{xsit))Sx{t) + cr-^fcii-^i (g) DG{xs{t - Tg))Sx{t - Tg) 

(59) 

+cr^fc2L^^ ® DG{xs{t - Tg))5x{t - Tg) + a^L^ ® DH{xs{t - Th))Sx{t - th). 

Following Sec. IIC, it can be shown that a necessary and sufficient condition for stability of the synchronous solution 
for the hypernetwork (|59p is that the maximum Lyapunov exponent of the low-dimensional equation, 

e^it) = DF{xS))dk{t) + a^ihXt - kia)DG{xs{t - Tg))6kit - Tg) + a'' \^ D H {x s{t - Th))ek{t - t^), (60) 

is negative for fc = 1, {N — 1), where a = J2f=i ^^'^ respectively the (complex) eigenvalues of the 

matrices and that are associated with the same eigenvectors, i.e., such that L^Vk = A^Wfc and L^Vk — Wfe. 
Note that the one eigenvalue A;^ — Xfj = is not relevant in determining stability. Now the question arises how 
likely it is that an arbitrary Laplacian matrix L"^ can be decomposed in the form (|58p . In general terms, an A^-squared 
matrix is determined by its N'^ entries. At the same time, we are allowed 2A^ degrees of freedom in the decomposition 
([55]). i.e., N degrees of freedom in choosing the entries ai, 02, ajv of the C-matrix L^^ and N degrees of freedom in 
choosing the eigenvalues of the matrix L"^^. It follows that only in the case that = 2, a decomposition in the form 
(1551) is guaranteed irrespective of the choice of the two Laplacian matrices and . This leads to the conclusion 
that the stability of the synchronous solution for an arbitrary N = 2-hypernetwork can always be associated with the 
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MLE of the low-dimensional equation (l60l) for k = 1. 
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